#include "TH3D.h"
#include "TRandom3.h"
#include "TVector3.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TF1.h"

const Double_t EPSILON = 1e-10;

void Deflection(TVector3 * vec, const Double_t dtheta, const Double_t dphi)
{
  /*
//test

.L dotest.C+
TH3D *hh = new TH3D("hh","",100,-5,5,100,-5,5,100,-5,5)
hh->Draw()
TVector3 aa(1,1,1)
const TVector3 a0=aa; 
const Int_t nn=50;
hh->Fill(aa.X(),aa.Y(),aa.Z());
for(int ii=0; ii<(3*nn); ii++){ Deflection(&aa, (TMath::PiOver4()),(TMath::TwoPi())/nn*ii); hh->Fill(aa.X(),aa.Y(),aa.Z()); aa=a0;}
hh->SetMarkerStyle(20);
hh->SetMarkerSize(1);
hh->Draw();

  */
  const Double_t theta0 = vec->Theta();
  const Double_t phi0 = vec->Phi();

  vec->SetMagThetaPhi(vec->Mag(), dtheta, dphi);

  vec->RotateY(theta0);
  vec->RotateZ(phi0);
}

//reasonable input for Rutherfordian scattering opt =1 and =2
//the dip at ~0 is due to regularization, the dip disappears at regpar->0
//dotest(1,0.000001)
//dotest(2,0.000001)
//BinLog x may be better
void dotest(const Int_t topt=0, const Double_t regpar=-999)
{
  TCanvas * cc = new TCanvas("cc","",800,800);
  cc->Divide(2,2);

  TH3D * hstep = new TH3D("hstep","", 50, -2, 2, 50, -2, 2, 50, -2, 2);
  TH3D * hsum = new TH3D("hsum","", 50, -2, 2, 50, -2, 2, 50, -2, 2);
  TH1D * htheta = new TH1D("htheta","", 180, 0, TMath::Pi());
  TH1D * hx = new TH1D("hx","", 100, -2,2);

  TRandom3 ran(123);

  const Int_t nexp = 10000;
  for(Int_t iexp=0; iexp<nexp; iexp++){
    TVector3 walker;
    walker.SetXYZ(0,0,1);

    const Int_t nstep = 1000;//ran.Poisson(1000);
    for(Int_t istep=0; istep<nstep; istep++){
      const Double_t dphi = ran.Rndm()*TMath::TwoPi();
      const Double_t rr = ran.Rndm();

      Double_t dtheta = -999;

      if(topt==0){
        //1) isotropic
        //(1-cos t)/2 = rr
        dtheta = TMath::ACos(1-2*rr);
      }
      else if(topt==1){
        //2) Rutherfordian
        //eq 8.46 Bielajew A
        dtheta = TMath::ACos( 1-2*regpar*(1-rr)/(regpar+2*rr) );
      }
      else if(topt==2){
        //3) Rutherfordian, small angle
        const Double_t tmparg = 2*regpar*rr/(1-rr);
        if(tmparg<EPSILON){
          continue;
        }
        dtheta = TMath::Sqrt(tmparg);
        if(dtheta>TMath::Pi()){
          continue;
        }
      }
      else{
        printf("wrong topt %d\n", topt); exit(1);
      }

      //cf Bielajew chapter 7, eq. 7.13
      Deflection(&walker, dtheta, dphi);

      //-------------
      TVector3 vecStep;
      vecStep.SetMagThetaPhi(1, dtheta, dphi);
      
      hstep->Fill(vecStep.X(), vecStep.Y(), vecStep.Z());
      //----------

      //constant delta_rho is fine if only deflection is concerned
      const Double_t rho = walker.Mag();
      walker.SetMagThetaPhi(rho+1, walker.Theta(), walker.Phi());
    }

    const TVector3 vecDel = walker;//-v0;
    htheta->Fill(vecDel.Theta());
    const Double_t nor = nstep;
    hx->Fill(vecDel.X()/nor);
    hsum->Fill(walker.X()/nor, walker.Y()/nor, walker.Z()/nor);
    //walker.Print();
  }

  cc->cd(1);
  hstep->Draw();

  cc->cd(2);
  hsum->Draw();

  cc->cd(3);
  htheta->Draw();

  TString ftit;
  Double_t pars[10];
  if(topt==0){
    ftit="[0]*sin(x)";
    pars[0]=1;
  }
  else if(topt>=1){
    ftit="[0]*sin(x)/pow(1-cos(x)+[1], 2)";
    pars[0]=1;
    pars[1]=regpar;
  }
  TF1 * ff=new TF1("ff", ftit,0,5);
  ff->SetParameters(pars);
  htheta->Fit(ff);

  cc->cd(4);
  hx->Draw();
  
  if(topt==0){
    ftit="[0]*TMath::Gaus(x,[1],[2],1)";
    pars[0]=1;
    pars[1]=0;
    pars[2]=1;
  }
  else if(topt>=1){
    //ftit="[0]*TMath::CauchyDist(x,[1],[2])";
    ftit="[0]*TMath::Gaus(x,[1],[2],1)";
    pars[0]=1;
    pars[1]=0;
    pars[2]=1;
  }
  TF1 * fx=new TF1("fx", ftit,-5,5);
  fx->SetParameters(pars);
  hx->Fit(fx);
}
